{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "positive-florence",
   "metadata": {},
   "source": [
    "# Problem 14.03 Reversible and Isothermal Compression of Liquid Water"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecological-electric",
   "metadata": {},
   "source": [
    "A flow of 2000 kg/h liquid water at 25 °C and 1 bar is pumped to a pressure of 100 bar. The pump is \"cooled\", so the process is reversible and isothermal. What is the duty of the pump shaft, and the energy that must be removed from the water being compressed?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tropical-beatles",
   "metadata": {},
   "source": [
    "## Solution\n",
    "\n",
    "We can use the high-accuracy IAPWS-95 implementation of the properties of water to easily and extremely accurately calculate these values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "advisory-execution",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shaft power is 5504.05633851 W\n",
      "The cooling duty is 431.1770 W\n"
     ]
    }
   ],
   "source": [
    "from scipy.constants import bar, hour\n",
    "import numpy as np\n",
    "from thermo import FlashPureVLS, IAPWS95Liquid, IAPWS95Gas, iapws_constants, iapws_correlations\n",
    "from scipy.integrate import quad\n",
    "from chemicals import property_molar_to_mass\n",
    "import numpy as np\n",
    "\n",
    "T1 = T2 = 25 + 273.15\n",
    "P1 = 1*bar\n",
    "P2 = 100*bar\n",
    "\n",
    "\n",
    "liquid = IAPWS95Liquid(T=T1, P=P1, zs=[1])\n",
    "gas = IAPWS95Gas(T=T1, P=P1, zs=[1])\n",
    "flasher = FlashPureVLS(iapws_constants, iapws_correlations, gas, [liquid], [])\n",
    "\n",
    "mass_flow = 2000/hour\n",
    "mole_flow = property_molar_to_mass(mass_flow, MW=iapws_constants.MWs[0])\n",
    "\n",
    "entry = flasher.flash(T=T1, P=P1)\n",
    "leaving = flasher.flash(T=T2, P=P2)\n",
    "\n",
    "def to_int(P, flasher):\n",
    "    state = flasher.flash(T=T1, P=P)\n",
    "    return state.V()\n",
    "integral_result = quad(to_int, P1, P2, args=(flasher,))[0]\n",
    "shaft_duty = integral_result*mole_flow\n",
    "\n",
    "cooling_duty = shaft_duty - (leaving.H() - entry.H())*mole_flow\n",
    "\n",
    "print(f'The shaft power is {shaft_duty:.8f} W')\n",
    "print(f'The cooling duty is {cooling_duty:.4f} W')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abstract-invite",
   "metadata": {},
   "source": [
    "The above shows the numerical integral calculation. That is the correct formulation.\n",
    "\n",
    "However, it can be a little unintuitive. We can contrast this with another calculation - a series of tiny isentropic compression, then cooling steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "knowing-queen",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The shaft power is 5504.0608 W\n",
      "The cooling duty is 431.1815 W\n"
     ]
    }
   ],
   "source": [
    "cooling_duty = 0\n",
    "compressing_duty = 0\n",
    "increments = 30 # Number of increments\n",
    "dP = (P2 - P1)/increments\n",
    "\n",
    "old_state = entry\n",
    "for i in range(increments):\n",
    "    P = P1+(i+1)*dP\n",
    "    \n",
    "    # Compress another increment of pressure\n",
    "    new_compressed_state = flasher.flash(S=old_state.S(), P=P)\n",
    "    compressing_duty += (new_compressed_state.H() - old_state.H())*mole_flow\n",
    "    \n",
    "    # Cool back to T1 at new pressure\n",
    "    new_cooled_state = flasher.flash(T=T1, P=P)\n",
    "    cooling_duty += (new_compressed_state.H() - new_cooled_state.H())*mole_flow\n",
    "    \n",
    "    old_state = new_cooled_state\n",
    "\n",
    "print(f'The shaft power is {compressing_duty:.4f} W')\n",
    "print(f'The cooling duty is {cooling_duty:.4f} W')"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}